library(rstan)

cat("\nCCES 2012 PID\n\n")

# Set to the local directory of the replication data
setwd("~/Downloads/Replication_Data/")
n_cores <- 8
options(mc.cores = n_cores)

L2012D <- readRDS("Data/CCES-L2012D.rds")
L2012R <- readRDS("Data/CCES-L2012R.rds")

##### Stack data to feed to Stan

model_data_2012D <- list(n_items = length(unique(L2012D$j)),
                        n_respondents = length(unique(L2012D$k)),
                        n_obs = nrow(L2012D),
                        k = L2012D$k, j = L2012D$j, y = L2012D$y)

model_data_2012R <- list(n_items = length(unique(L2012R$j)),
                        n_respondents = length(unique(L2012R$k)),
                        n_obs = nrow(L2012R),
                        k = L2012R$k, j = L2012R$j, y = L2012R$y)

##### Compile to C

model_am <- stan_model(file = "AM_Sigma.stan")

##### Sample from posterior

exclude_pars <- c("beta_raw", "alpha_raw", "zeta_raw")

# Democrats
n_items <- length(unique(L2012D$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2012D <- sampling(model_am, data = model_data_2012D,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2012D, "Fitted_Models/posterior_cces_2012D.rds", compress = "xz")


# Republicans
n_items <- length(unique(L2012R$j))
inits <- replicate(n_cores, list(list(zeta_raw = c(-1, 1, -1, 1, -1, runif(n_items-7, -2, 2)))))

posterior_cces_2012R <- sampling(model_am, data = model_data_2012R,
                                warmup = 1000, iter = 1250, seed = 1,
                                include = FALSE, pars = exclude_pars,
                                init = inits,
                                refresh = 50, save_warmup = FALSE,
                                control = list(adapt_delta = 0.8,
                                               max_treedepth = 10,
                                               stepsize = 1),
                                open_progress = TRUE, thin = 1, chains = n_cores)

saveRDS(posterior_cces_2012R, "Fitted_Models/posterior_cces_2012R.rds", compress = "xz")

